rm(list=ls())

library(cherry)
library(hommel)
library(globaltest)
library(ctgt)


rejectedpathway <- function(data, pathinfo){
  y = data[,ncol(data)]
  X = data[,-ncol(data)]
  ##globaltest
  mygt <- function (hps){
    sqrW = sqrt(mean(y)*(1-mean(y)) ) # sqrt of covariance of y
    WIHZ = sqrW *(sweep(X,2,colMeans(X))) ## W^{1/2}*(I-H)Z
    IHZ = WIHZ/sqrW ##(I-H)Z
    # the full model
    test = sum(colSums(y*IHZ[,hps,drop=F])^2)# full test statistic = y^t (I-H) Z Z^t (I-H) y
    lamt = round(eigen(tcrossprod(WIHZ[,hps,drop=F]),symmetric = T,only.values = T)$values,8)
    pv(test,lamt)
  }
  
  ## DAG
  stc = construct(pathinfo)
  DAG <- DAGmethod(stc, mygt,isadjusted = T)
  resdag  = pathinfo[!is.na(pvalue(DAG))]
  
  ## SH
  SH = tryCatch(structuredHolm(stc, mygt,isadjusted = T), error = function(e) 0) 
  if(class(SH)!="DAG"){
    sher = rep(FALSE, length(pathinfo))
  } else{
    sher = !is.infinite(pvalue(SH)) 
  }
  ressh = pathinfo[sher]
  
  ## FL
  names(pathinfo) = pathinfo
  fc = findFocus(pathinfo, maxsize = 10, atoms = TRUE)
  resf = focusLevel(mygt, pathinfo, focus = fc)
  pf = unlist(resf$focuslevel)
  resfl = pathinfo[pf<0.05]
  
  ## Hommel
  xs = colnames(X)
  pval = sapply(xs, mygt )
  hom=hommel(pval,simes = T)
  pa = sapply(pathinfo, function(x) localtest(hom, x,tdp=0))
  reshommel = pathinfo[which(pa<0.05)]
  
  ##ctgt
  resctgt20t = sapply(pathinfo, function(i) actgt(y,X,xs, as.character(i),maxit= 20000)[1])
  
  res = list(ctgt = pathinfo[resctgt20t=="reject"], dag = resdag, sh = ressh,fl = resfl, ctst = reshommel )
  return(res)
}

#####################################################################
### Eisner
#####################################################################
load("Eisner/d_uri.RData")
load("Eisner/allpathwayuri.RData")
res_eisner = rejectedpathway(data = d_uri, pathinfo = allpathwayuri)

save(res_eisner, file="Eisner/res_eisner.RData")

#####################################################################
### Bordbar
#####################################################################
load("Bordbar/d23.RData")
load("Bordbar/allpathway23.RData")
res_bordbar = rejectedpathway(data = d23, pathinfo = allpathway23)

save(res_bordbar, file="Bordbar/res_bordbar.RData")

#####################################################################
### Taware
#####################################################################
load("Taware/d760.RData")
load("Taware/allpathway760.RData")
res_taware = rejectedpathway(data = d760, pathinfo = allpathway760)

save(res_taware, file="Taware/res_taware.RData")

#####################################################################
### Al-Mutawa
#####################################################################
load("Al-Mutawa/d541.RData")
load("Al-Mutawa/allpathway541.RData")
res_almutawa = rejectedpathway(data = d541, pathinfo = allpathway541)

save(res_almutawa, file="Al-Mutawa/res_almutawa.RData")
